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Abstract 

We derive a master stability function (MSF) for synchronization in networks of coupled dynami- 
cal systems with small but arbitrary parametric variations. Analogous to the MSF for identical sys- 
tems, our generalized MSF simultaneously solves the linear stability problem for near-synchronous 
states (NSS) for all possible connectivity structures. We also derive a general sufficient condi- 
tion for stable near-synchronization and show that the synchronization error scales linearly with 
the magnitude of parameter variations. Our analysis underlines significant roles played by the 
Laplacian eigenvectors in the study of network synchronization of near-identical systems. 
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Introduction — Synchronization in its various forms has been a highly popular and excit- 
ing developing topic in the recent literature on chaotic oscillators [HQS]. Applications have 
ranged widely from biology [U[20] to mathematical epidemiology [18], and modelling animal 
gaits jH] to engineering of communications devices pfl |3] , including many developments in 
complex networks (see e.g., Refs. [TUJ [2H 1221 [221 12S] and a review [IS])- However, the pre- 
ponderance of the work has focused on identical synchronization since it is in this situation 
whereby a complete analysis can be carried forward by the master stability formalism de- 
veloped in the seminal work [12] . While other forms of synchronization are discussed in the 
literature, of particular interest here is nearly-synchronous state behavior of the systems that 
are slightly detuned from identical synchronization, which may or may not be associated 
with an invariant manifold [13] normally required to describe generalized synchrony El E] ■ 

In this letter we consider a coupled dynamical system consisting of N units coupled 
through some underlying network. The equations of motion reads: 

N 

Wi = f(wi,fii) - g^LijHiwj), i = 1,2, . . . ,JV, (1) 

3=1 

where / : 9ft mxp — > 9ft m is the parameterized dynamics of an isolated unit; Wi G 9ft m is the 
dynamical variable for the ith unit; /ij G 3ft p is the corresponding parameter; L G dt NxN is 
the graph Laplacian [2S]; H : 3ft m — > W a is a uniform coupling function; and g G 9ft is the 
uniform coupling strength (usually > for diffusive coupling). 

Note that we can represent the whole system conveniently by using Kronecker product 
representation: 

w = f(w,fi)-g-L®H(w), (2) 

where w= [w^, ...wJj] T is a column vector of all the dynamic variables, and likewise for 
fj, and /; and ® is the usual Kronecker product (or direct product) [2]. 

System ([T| has been studied mostly in the case in which the parameter /Zj is the 
same for each individual oscillator, often resulting in identical synchronization where 
maxjj ||ii?i(£) — Wj(t)\\ — ► as t — ► oo. The stability of such states can be analyzed by 
master stability functions (MSF) |12j . 

However, a noiseless system with exactly the same parameters is impossible in practice. 
It is known that parameter mismatch among the individual oscillators can cause bursts due 
to the instability of typical periodic orbits embedded in the synchronized chaotic attractor 
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[2T] ; even within a stable region where no bubbling will occur, the states of different units 
will still not approach exactly the same function of time, but instead come close to each 
other within a neighborhood of the identical synchronization state [21]. This phenomena 
was first reported in [1] for two coupled Lorenz oscillators, where the variations of individual 
units from the identical synchronization manifold was found to scale linearly with respect 
to the magnitude of parameter mismatch when the mismatch is small. In [2T], a variational 
equation analogous to our Eq. ^ was used to study the progressive loss of synchronization 
stability due to bursting, which is also a relevant and interesting phenomenon. In this 
letter we develop an extended master stability framework for systems with near-identical 
parameters and derive stability conditions for stable near-synchronization. 

Near-Synchronous State (NSS) — Assume that the parameters pi in Eq. (JTJ are close to 
each other and do not change with time. Let the average parameter be p = 4 J2f=i Pi an d 
the parameter mismatch be 8 fa = fa — p. With appropriate choices of coupling strength 
g and network structure L, the system can have a near synchronous state (NSS) in which 
maxjj ||u>i(t) — iWj* (f ) 1 1 < c as t — > oo for some small constant c > 0. When the system un- 
dergoes such near-synchronization, the trajectories of individual units are well approximated 
by the average trajectory w = 4 J2 i=1 w i7 which is governed by 



where we have defined dj = ~^J2i=i^ij [211 • With this equation, we can discuss dynamics 
of the bulk, or coarse scale behavior. 

Inhomogeneity in Variational Equations — Define the variation on each individual unit 
to be rji = Wi — w for % = 1, 2, N. The variational equations is then 



Assuming that the variations rji and the parameter mismatch S fa are small, we expand 
around w and p to obtain 




(3) 




N 




(4) 



N 




+D l ,f(w,fi)Sfa. 



(5) 
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We have used Y^f=i Vj = Y,jLi Wj - N ■ w — 0, Y!,jLi = J2j=i A*j — N ■ ft — 0, and 

>7V 



J^j=i 4? = Si j Lij = in the derivation. Putting all the and 5 fa in column vectors r] and 
8 ft, respectively, and omitting the arguments (w,fl) for simplicity, we obtain the variational 
equation for the NSS: 



IN 



D w f — g ■ G ® DH t) + I N ®D„f 5fi, (6) 



where the modified graph Laplacian G is defined by G = L — [1, 1, l] r • [cfi, <i 2 , ...djv] and 
In is the N x N identity matrix. 

Interestingly, the eigenvalues of G, Ai, Ag, • • • , Xn are exactly the same as those of L, 
and the vector [1,1,..., 1] is the eigenvector of both L and G associated with Ai = 0. 
Furthermore, any other eigenvector v' of G associated with eigenvalue A can be obtained 
from the corresponding eigenvector v of L by the transformation v ' = v — [v, ...,v] T , which 
simply shifts each component of v by a constant v = jYlf=i^j v j- More importantly, if 
there exist diagonalization forms L = QAQ^ 1 and G = PAP" 1 , then the corresponding 
rows of Q -1 and P -1 (the left eigenvectors of L and G, respectively) are parallel to each 
other, except for the first rows that correspond to Ai = 0. 

When all the parameters fa are the same, the second term in the Eq. ^ disappears, and 
what is left is a homogeneous ODE system for 77, which may be diagonalized to obtain an 
equation analogous to the well-known master stability equation [12], with the only difference 
that here we have a modified graph Laplacian G. Interestingly, in the case of no parameter 
mismatch, this difference would not lead to different conclusions since the stability analysis 
depends on the graph structure only through the Laplacian eigenvalues, not eigenvectors. 

We now focus on the case in which, if there were no parameter mismatch, the system 
would undergo stable identical synchronization, i.e., the variation 77 would go to zero asymp- 
totically. This situation occurs if the system represented by /, H, L and g are in the stable 
regime [12]. Because of the inhomogeneous part [pv ® D^f^dfi due to parameter mismatch, 
the variational system ^ in general may not be asymptotically stable. We will show, how- 
ever, that when the parameter mismatch is small, there may exist a NSS where r\ stays close 
(although not equalling) to zero. Indeed, we will show that the variational system is stable 
(i.e. the solution 77 is bounded as t — > 00) and the bound for the solution depends linearly 
on the norm of the parameter mismatch 8 ft. 

Extended Master Stability Equation and Function — We may uncouple the variational 



equation by diagonalizing the modified graph Laplacian G: G = PAP 1 [28] for some 
invertible matrix P. Making the change of variable £ = (P _1 ® I m )r], we obtain 



c 



In ® D w f - g- k®DH C+ P~ l ® DJ 8fi. 



(7) 



The homogeneous part in Eq. ^ has block diagonal structure and we may write for each 
eigenmode i > 2 

N 



D w f - gXiDH Q + D^f ■ w^tf/i 



(8) 



where ity is the jth component of the ith row in the matrix P , i.e., Ui is the zth left 



eigenvector of G. The vector YljLi u ij^^j 1S the weighted average of parameter mismatch 
vectors, with the weights given by the components of the left eigenvector associated with 
Aj. It may also be thought of as an inner product of the parameter mismatch vector and 
the corresponding left eigenvector. We comment here that if one uses the original graph 
Laplacian L instead, the resulting equation would be equivalent to Eq. (|8j, since the spectra 
of L and G are the same and corresponding left eigenvectors are parallel except for those 
associated with Ai = 0. 

From Eq. ([8]), we define an extended master stability equation for near identical coupled 
dynamical systems: 



D w f — a ■ DH £ + D„f ■ ip 



(9) 



where we have introduced two auxiliary parameters, a (complex) scalar a and ip G W. 
Once the stability of Eq. ^ is determined as a function of a and ip, the stability of the ith 
eigenmode can be found by simply setting a = g\i and ip = Uij^j- The problem is 

thus decomposed into two separate parts: one that depends only on the individual dynamics 
and the coupling function, and the other that depends only on the graph Laplacian and 
parameter mismatch. Note that the latter not only depends on the spectrum of L as in 
|12j . but also on the combination of the left eigenvectors and parameter mismatch. Thus, 
we have reduced the stability analysis of the original miV-dimensional problem to that of 
m-dimensional problem with one additional parameter, combined with an eigenproblem. 

Note that to analyze the stability of the original system using the master stability equa- 
tion, we need the associated average trajectory w, which can only be obtained by solving the 
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original system, and is impractical for large networks. We found, however, that in practice 
(as we will confirm in an example below) one may instead use a trajectory s of a single 
auxiliary average unit: s = f(s,p,). We conjecture that under suitable conditions on the 
system, the trajectory s of the average unit shadows the average trajectory w [29|. 

The associated master stability function Q(a,ip) is then defined to be the asymptotic 
value of the norm of £ as a function of a and ip, given that a leads to asymptotic stable 
solution of the homogeneous part. In the case of symmetrically coupled networks, for which 
G = L is symmetric, the matrix P can be chosen to be orthogonal, allowing us to predict 
the square-sum synchronization error in the original system ([I]) from Q(a, tp): 



N N N 



i=l i=2 i=2 

where and ipi correspond to the iih eigenmode and 1 1 ■ 1 1 denotes the Euclidean norm. 

Conditions for Stable Synchronization — In the previous section we have derived a generic 
stability equation ^ for analyzing the stability of synchronization of coupled dynamical 
system ([!]) . To analyze the stability, we now assume that the largest Lyapunov exponent of 
the synchronous trajectory associated with the homogeneous variational equation 



D w f - aDH 



£ (11) 



is negative for a given a, so that without parameter mismatch the error mode corresponding 



to this specific a goes to zero exponentially. In this case, the solution of Eq. (11) can 
be written as: £*(£) = $(£, 0)£(0), where $(t, r) is the fundamental transition matrix [30J, 
satisfying 

\\®(t,T)\\ < ie - x{t - T) (12) 

for t > t and some finite positive constants 7 and A. We should note that in the case of 
generalized synchrony, the loss of stability of the invariant manifold need not proceed mono- 
tonically and uniformly in space. It is known that parameter mismatch can cause bursting 
due to increasing instability of embedded transversely unstable periodic orbits which cause 
short-time positivity of Lyapunov exponents [TTJ, &T\, and this can be correspondingly inter- 
preted from Eq. (12). Such transition has been called bubbling bifurcation [TUl [TTJ due to 
basin riddling. 

The solution to Eq. ^ can then be expressed by 

m = Ht, 0)£(0) + / $(t, r)b(r)dr } (13) 
Jo 



where 6(r) = D^f^sij), /x) • ■0. Under the condition of Eq. (12), we can show that given 



by Eq. ( 13 ) is bounded by the following inequality 

wmw < \m,o)\\ -iie(o)n + 



< 7e 
7 



-At 



Hewn 



|$(t,r)||dr-sup||6(t)|| 
o t 



sup ||6(0|| 



i 



sup ||6(t)|| as £ — > oo. 
A t 



(14) 



Thus, the inhomogeneous master stability equation is stable, i.e., the solution to Eq. (j9|) is 
bounded asymptotically as long as i) the homogeneous system is exponentially stable, or 
equivalently, the maximal Lyapunov exponent is negative; and ii) the inhomogeneous part 
6(t) = D^f{s{r),Jj) ■ ip is bounded. 



Eq. (12) and Eq. (13) also allow us to analyze quantitatively the magnitude of asymptotic 



error of a near-identical system. If the magnitude of parameter mismatch is scaled by a factor 



c, keeping all other parameters fixed, it follows from Eq. (13) that the corresponding solution 
will be 

£(£) = $(£,0)£(0)+c [ $(t,T)b(r)dT, 



(15) 
Jo 

where £(£) denotes the variation evolution of the original unsealed near-identical system. 



Now the first term of both Eq. (13) and Eq. (15) goes to zero exponentially according to 



Eq. (12), so that asymptotically we have = c£(£), i.e., the variation is scaled by the 
same factor correspondingly. 

Examples of Application — We consider each individual unit w = [x,y, z] T governed by 
the Lorenz equations: 

x = a(y-x), 
y = x (r - z) - y, 

z = xy — (3z, (16) 

where parameters a — 10, /3 — |, and we consider mismatch between units in r, i.e., r 
corresponds to /i in Eq. Q. So we have 



-a a 

r — z —1 —x 
V x -(3 



(17) 



7 



45 

a 



-0.04 



'0 



FIG. 1: Density plot of extended master stability function Q(a, ip) associated with arbitrary 



networks of near-identical Lorenz systems. It is estimated by y y f \\£(t)\\ 2 dt with T = 200 (||.|| 
denotes the Euclidean norm), where is obtained by numerically integrating Eq. (j9j with a time 
step of 0.001 and discarding initial transient. Here we have used the coupling function H{x) = x. 

and D^f = [0, x, 0] T . The coupling function H is taken to be H(x) = x, so that DH(s) = J 3 
(Vs). With these choices of / and H, we numerically integrate Eq. ^ for a range of a 
and ij) and estimate the asymptotic norm of which gives Q(a, ip) shown in Fig. [TJ As 
shown in Fig. [2] for a 4-node network example in the inset, this estimated ip), combined 
with Eq. ( JToj ) gives fairly good predictions for the actual synchronization error in the full 
system 0. In addition, Fig. [2] confirms that the actual synchronization error scales linearly 
with the magnitude of the parameter mismatch, as predicted by our analysis. 

Summary and Discussion — In this letter we have analyzed the stability of synchroniza- 
tion in a network of coupled near-identical dynamical systems. We have shown that the 
well-known master stability approach can be extended to this general case, allowing us to 
solve the part of the problem that depends on the individual node dynamics, independently 
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0.02 

FIG. 2: Comparison of predicted synchronization error with actual error. For the 4-node network 
shown in the inset, the error prediction Yli=2 ^( a «; ''Pi) 2 (dashed line) was computed using Q dis- 
played in Fig. [TJ The values of a« and ip{ (shown in Fig. [l]by arrows for e = 0.01) were obtained 
from the Laplacian eigenstructure and the parameter mismatch pattern indicated in the inset as 
a function of e. Actual error (squares) was estimated by ^ Jq ^i=i H^WII 2 ^ with T = 200 com- 
puted from numerical integration of the full system ([!]) after discarding initial transient. We used 
g = 5 in all calculations. 

of the network structure and the parameter mismatch pattern over the network. We have 
demonstrated the validity of our analysis using a small example network of coupled Lorenz 
systems. The extended MSF gives simplified, accurate, and practical estimate of the magni- 
tude of variation in a near-identical system, provided that the corresponding identical system 
undergoes stable synchronization according to the original MSF analysis. Furthermore, our 
results highlight the relevance of the Laplacian eigenvector structure, in addition to the full 
eigenvalue spectrum, in determining the amount of dynamical variation due to parameter 
mismatch among individual dynamics. This suggests that detailed knowledge of the graph 
structure may be important for the design of robust and reliable systems. 
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